A short description of the post.
packages = c('maptools', 'sf', 'raster', 'spatstat', 'tmap', 'tidyverse', 'plotly', 'ggthemes')
for (p in packages){
if(!require(p, character.only = T)){
install.packages(p)
}
library(p,character.only = T)
}
Importing Shapefile using st_read() of sf package.The output object is in tibble sf object class.
mpsz_sf <- st_read(dsn = "data/shapefile", layer = "MP14_SUBZONE_WEB_PL")
Reading layer `MP14_SUBZONE_WEB_PL' from data source
`C:\darrylkwok\IS415_blog\_posts\2021-09-13-in-class-exercise-5\data\shapefile'
using driver `ESRI Shapefile'
Simple feature collection with 323 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
Projected CRS: SVY21
mpsz_sf
Simple feature collection with 323 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
Projected CRS: SVY21
First 10 features:
OBJECTID SUBZONE_NO SUBZONE_N SUBZONE_C CA_IND
1 1 1 MARINA SOUTH MSSZ01 Y
2 2 1 PEARL'S HILL OTSZ01 Y
3 3 3 BOAT QUAY SRSZ03 Y
4 4 8 HENDERSON HILL BMSZ08 N
5 5 3 REDHILL BMSZ03 N
6 6 7 ALEXANDRA HILL BMSZ07 N
7 7 9 BUKIT HO SWEE BMSZ09 N
8 8 2 CLARKE QUAY SRSZ02 Y
9 9 13 PASIR PANJANG 1 QTSZ13 N
10 10 7 QUEENSWAY QTSZ07 N
PLN_AREA_N PLN_AREA_C REGION_N REGION_C
1 MARINA SOUTH MS CENTRAL REGION CR
2 OUTRAM OT CENTRAL REGION CR
3 SINGAPORE RIVER SR CENTRAL REGION CR
4 BUKIT MERAH BM CENTRAL REGION CR
5 BUKIT MERAH BM CENTRAL REGION CR
6 BUKIT MERAH BM CENTRAL REGION CR
7 BUKIT MERAH BM CENTRAL REGION CR
8 SINGAPORE RIVER SR CENTRAL REGION CR
9 QUEENSTOWN QT CENTRAL REGION CR
10 QUEENSTOWN QT CENTRAL REGION CR
INC_CRC FMEL_UPD_D X_ADDR Y_ADDR SHAPE_Leng
1 5ED7EB253F99252E 2014-12-05 31595.84 29220.19 5267.381
2 8C7149B9EB32EEFC 2014-12-05 28679.06 29782.05 3506.107
3 C35FEFF02B13E0E5 2014-12-05 29654.96 29974.66 1740.926
4 3775D82C5DDBEFBD 2014-12-05 26782.83 29933.77 3313.625
5 85D9ABEF0A40678F 2014-12-05 26201.96 30005.70 2825.594
6 9D286521EF5E3B59 2014-12-05 25358.82 29991.38 4428.913
7 7839A8577144EFE2 2014-12-05 27680.06 30230.86 3275.312
8 48661DC0FBA09F7A 2014-12-05 29253.21 30222.86 2208.619
9 1F721290C421BFAB 2014-12-05 22077.34 29893.78 6571.323
10 3580D2AFFBEE914C 2014-12-05 24168.31 30104.18 3454.239
SHAPE_Area geometry
1 1630379.3 MULTIPOLYGON (((31495.56 30...
2 559816.2 MULTIPOLYGON (((29092.28 30...
3 160807.5 MULTIPOLYGON (((29932.33 29...
4 595428.9 MULTIPOLYGON (((27131.28 30...
5 387429.4 MULTIPOLYGON (((26451.03 30...
6 1030378.8 MULTIPOLYGON (((25899.7 297...
7 551732.0 MULTIPOLYGON (((27746.95 30...
8 290184.7 MULTIPOLYGON (((29351.26 29...
9 1084792.3 MULTIPOLYGON (((20996.49 30...
10 631644.3 MULTIPOLYGON (((24472.11 29...
Projection is in SVY21.
read_rds() of readr package is used instead of readRDS() of base R. This is because the output of read_rds() is a tibble object.
CHAS <- read_rds("data/rds/CHAS.rds")
childcare <- read_rds("data/rds/childcare.rds")
Note that there are some issues found in the childcare dataframe, because the Lat and Lng attributes should be in numeric data type. The coordinate fields seems to be in decimal degree. Therefore, this results in having assumed the WGS84 referencing system.
st_crs(CHAS)
Coordinate Reference System: NA
CHAS_sf <- st_as_sf(CHAS, coords = c("X_COORDINATE", "Y_COORDINATE"), crs = 3414)
Note: st_as_sf accepts coordinates that are in character data type, and converts it into the appropriate data type.
We see that it is in decimal degree in the childcare data, therefore we know that it is 4326 projection system and we will transfrom it to 3414
childcare_sf <- st_as_sf(childcare, coords = c("Lng", "Lat"), crs = 4326) %>%
st_transform(crs=3414)
This is how we can check the projection system. If the points are plotted on somewhere random, it means that the projection system is wrong.
All the darker points show that you have more than 1 childcare centers
tmap_mode("view")
tm_shape(childcare_sf) +
tm_dots(alpha = 0.4, col = "blue", size = 0.05) +
tm_shape(CHAS_sf) +
tm_dots(alpha = 0.4, col = "red", size = 0.05)
childcare <- as_Spatial(childcare_sf)
CHAS <- as_Spatial(CHAS_sf)
mpsz <- as_Spatial(mpsz_sf)
Using *as.ppp() of maptools** package
childcare_sp <- as(childcare, "SpatialPoints")
CHAS_sp <- as(CHAS, "SpatialPoints")
mpsz_sp <- as(mpsz, "SpatialPolygons")
Using as.ppp() of maptools package
childcare_ppp <- as(childcare_sp, "ppp")
CHAS_ppp <- as(CHAS_sp, "ppp")
If you check childcare_ppp and CHAS_ppp for duplicates using the any() and duplicated() function, you will see that it will return TRUE
childcare_ppp_jit <- rjitter(childcare_ppp, retry=TRUE, nsim=1, drop=TRUE)
any(duplicated(childcare_ppp_jit))
[1] FALSE
CHAS_ppp_jit <- rjitter(CHAS_ppp, retry=TRUE, nsim=1, drop=TRUE)
any(duplicated(CHAS_ppp_jit))
[1] FALSE
Need a comma at the end to select all the columns since the statement before the comma is to select the rows
pg <- mpsz[mpsz@data$PLN_AREA_N == "PUNGGOL",]
pg_sp <- as(pg, "SpatialPolygons")
pg_owin <- as(pg_sp, "owin")
Extract out the childcares and clinics that are within punggol.
If you plot childcare_ppp_jit, you will only get the points. After extracting the points using the owin object, you will get the Punggol’s Boundary and the points within Punggol
childcare_pg <- childcare_ppp_jit[pg_owin]
CHAS_pg <- CHAS_ppp_jit[pg_owin]
plot(childcare_pg)

L_childcare <- envelope(childcare_pg, Lest, nsim=99, rank=1, global=TRUE)
Generating 99 simulations of CSR ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70,
71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99.
Done.
title <- "Pairwise Distance: L function"
Lcsr_df <- as.data.frame(L_childcare)
colour=c("#0D657D","#ee770d","#D3D3D3")
csr_plot <- ggplot(Lcsr_df, aes(r, obs-r))+
# plot observed value
geom_line(colour=c("#4d4d4d"))+
geom_line(aes(r,theo-r), colour="red", linetype = "dashed")+
# plot simulation envelopes
geom_ribbon(aes(ymin=lo-r,ymax=hi-r),alpha=0.1, colour=c("#91bfdb")) +
xlab("Distance r (m)") +
ylab("L(r)-r") +
geom_rug(data=Lcsr_df[Lcsr_df$obs > Lcsr_df$hi,], sides="b", colour=colour[1]) +
geom_rug(data=Lcsr_df[Lcsr_df$obs < Lcsr_df$lo,], sides="b", colour=colour[2]) +
geom_rug(data=Lcsr_df[Lcsr_df$obs >= Lcsr_df$lo & Lcsr_df$obs <= Lcsr_df$hi,], sides="b", color=colour[3]) +
theme_tufte()+
ggtitle(title)
text1<-"Significant clustering"
text2<-"Significant segregation"
text3<-"Not significant clustering/segregation"
# the below conditional statement is required to ensure that the labels (text1/2/3) are assigned to the correct traces
if (nrow(Lcsr_df[Lcsr_df$obs > Lcsr_df$hi,])==0){
if (nrow(Lcsr_df[Lcsr_df$obs < Lcsr_df$lo,])==0){
ggplotly(csr_plot, dynamicTicks=T) %>%
style(text = text3, traces = 4) %>%
rangeslider()
}else if (nrow(Lcsr_df[Lcsr_df$obs >= Lcsr_df$lo & Lcsr_df$obs <= Lcsr_df$hi,])==0){
ggplotly(csr_plot, dynamicTicks=T) %>%
style(text = text2, traces = 4) %>%
rangeslider()
}else {
ggplotly(csr_plot, dynamicTicks=T) %>%
style(text = text2, traces = 4) %>%
style(text = text3, traces = 5) %>%
rangeslider()
}
} else if (nrow(Lcsr_df[Lcsr_df$obs < Lcsr_df$lo,])==0){
if (nrow(Lcsr_df[Lcsr_df$obs >= Lcsr_df$lo & Lcsr_df$obs <= Lcsr_df$hi,])==0){
ggplotly(csr_plot, dynamicTicks=T) %>%
style(text = text1, traces = 4) %>%
rangeslider()
} else{
ggplotly(csr_plot, dynamicTicks=T) %>%
style(text = text1, traces = 4) %>%
style(text = text3, traces = 5) %>%
rangeslider()
}
} else{
ggplotly(csr_plot, dynamicTicks=T) %>%
style(text = text1, traces = 4) %>%
style(text = text2, traces = 5) %>%
style(text = text3, traces = 6) %>%
rangeslider()
}